"""
Module containing function to calculate the one-dimensional
psuedo steady-state drying model.
Created: Brent Maranzano
Last Modified: May 22, 2014
"""

import numpy as np
from scipy.optimize import fmin
from matplotlib import pyplot as plt
from pdb import set_trace

def vapor_pressure(T, A, B, C):
    """
    Calculate the vapor pressure at temperature, T, using
    parameters for the Antoine eq. and return the vapor pressure in Pa.
    T - temperature to calculate the vaport pressure (C)
    parameters - Antoine parameters 
    return P_vap - vapor pressure at T (Pa) 
    """
    P_vap = 10**(A-B/(C+T))
    P_vap = P_vap * 1.01325*10**5/760
    return P_vap

def mass_flux(z, P, z_c=0.065, P_v=7000., K=0.8*10**-12,
        h=8.6*10**-6, r=0.068):
    """
    Mass evaporation rate per unit cross section area, which is 
    equal equal to the mass flux rate given by Darcy's Law.
    velocity = -(permeability/viscosity)*grad(pressure)
    P : pressure at z (Pa)
    z : axial distance from the cake bottome (m)
    P_v : vessel pressure (Pa)
    z_c : cake height (m)
    K : vapor phase permeability for Darcy's law (m^2)
    h : vapor phase viscosity (Pa s)
    r : vapor phase density (kg/m^3)
    return j : mass flux (kg / m^2 s)
    """
    j = - r * (K/h) * (P_v - P) / (z_c - z)
    return j

def Newton_cooling_parameter(X, h_w=10., h_f=10., X_cr1=0.1,
        X_cr2=0.08):
    """
    Calculate Newton's cooling parameter, h.
    X : actual solvent concentration
    h_w : heat transfer coefficient from wall to cake (W / m^2 /K)
    h_f : heat transfer coefficient from heating fluid to wall (W / m^2 /K)
    X_cr1 : first critical solvent concentration
    X_cr2 : second critical solvent concentration
    return h : Newton cooling parameter (W / m^2 /K)
    """
    if X < X_cr1 and X > X_cr2:
        h_g = 1 / (1/h_f + 1/h_w)
        h = (h_f - h_g)/(X_cr1 - X_cr2) * (X - X_cr2) + h_g
    else:
        h = h_w
    return h

def jacket_temperature(t, Tf_i=20, grad_Tf=0.1, Tf_s=45.):
    """
    Calculate the jacket temperature for transient heat changer fluid
    temperature. #TODO Where did this calculation come from?
    t : time (s)
    Tf_i : initial heat exchanger fluid temperture (K)
    grad_Tf : slope of change of heat fluid per time (K / s)
    return Tf_s : set point heat exchanger temperature (K)
    """
    T_f = numpy.min([Tf_i + grad_Tf*t,Tf_s])
    return T_f

def heat_flux(T, z, h, T_f=293., l_d=0.169):
    """
    Calculation of the energy flux through the wall into the
    filter cake. 
    T : temperature of front (K)
    z : location of front (m) # TODO check this
    T_f : temperature of fluid (K)
    h : Newton cooling parameter (W / m^2 /K)
    l_d : thermal conductivity dry layer (W / m /K)
    return q_w : heat flux (W / m^2)
    """
    q_w = (1 / (1/h + z/l_d)) * (T_f - T)
    return q_w

def heat_capacity(X, cp_s=837., r_s=2100., cp_l=4186., r_l=985.):
    """
    Calculate the heat capacity of the solid solvent mixture
    per unit area of vessel cross section
    X - mass fraction of solvent
    cp_s - heat capacity solid (J/kg)
    r_s - density solid (kg/m^3)
    cp_l - heat capacity liquid (J/kg)
    r_l - density liquid (kg/m^3)
    return heat_capacity : heat capacity of solid/solvent mixture (J/kg)
    """
    heat_capacity = cp_l*r_l*(X/(X+1)) + cp_s*r_s*(1/(X+1))
    return heat_capacity

def local_temp_objective_function(T, params):
    """
    Objective funtion to be minimized to solve for local temperature.
    heat_flux(T) = mass_flux(P_sat(T))*H_vap + rate change_sensible_heat(T_i-T)
    [(mass_flux*H_vap + change_sensible_heat(T_i-T) - heat_flux)(T)]^2 = 0
    T - temperature to solve for steady state (K)
    params : dictioary of parameters used in the simulation
    """
    P_sat = vapor_pressure(T, A, B, C)
    j = mass_flux(P_sat, z, z_c, P_v, K, v, r)
    q = heat_flux(T, z, h, T_f, l_d)
    S = heat_capacity(cp_s, r_s, cp_l, r_l, X) * (z_c-z) * (T-T_i) / dt
    objective_function = (j * H_vap + S - q)**2
    return objective_function

def regression_diagnostics_inspection(self,T_array,P_parameters,z,P_v,K,v,r,h,T_f,H_vap,z_c,l_d):
    """
    function to return the set of variables: T, P_sat, j*H_vap, q,
    objective function to inspect the nonlinear regression calculation
    """
    data = numpy.zeros(T_array.shape[0],dtype=[('T','f4'),('P_sat','f8'),('j','f8'),
           ('q','f8'),('m','f8'),('objective','f8')])
    for i in range(T_array.shape[0]):
        T = T_array[i]
        P_sat = self.vapor_pressure(T,P_parameters)
        j = self.mass_flux(P_sat,z,z_c,P_v,K,v,r)
        q = self.heat_flux(T,z,h,T_f,l_d)
        objective_function = (j*H_vap-q)**2
        data['T'][i] = T
        data['P_sat'][i] = P_sat
        data['j'][i] = j
        data['q'][i] = q
        data['m'][i] = j*H_vap
        data['objective'][i] = objective_function
    return data

def simulate(r_s=2100., Y=0.45, l_d=0.169, cp_s=837.,
        r_l=985., cp_l=4186.,
        v=8.6*10**-6, mw=30., H_vap=9.186*10**5,
        p1=8.20417, p2=1642.89, p3=230.3,
        K=0.8*10**-12, h_w=10, h_f=10.,
        X_cr1=0.1, X_cs2=0.08,
        z_c=0.0065, D=0.003, P_v=7000., Tf_i=293., grad_Tf=0.0002, Tf_s=318.,
        To=293., Xo=0.2, dt=20.): 
    """
    Psuedo steady state Agitated Filter Drying model in one dimension
    See Kohout et al Chem. Eng. Sci. 61 (2006)
    Developed by: Brent Maranzano
    Created date: 18, May 2011
    Last modified date: 27, June 2011
        r_s : solid density (kg/m^3)
        Y : dry solid porosity 
        l_d : thermal conductivity dry layer (W/m/K)
        cp_s : solid heat capacity (J/kg K)
        r_l : liquid density (kg / m^3)
        cp_l : liquid heat capacity (J / kg K)
        v : vapor viscosity (Pa s)
        mw : vapor molecular weight (g/mol) 
        H_vap : enthalpy of vaprization (J/mol) #TODO check units
        p1, p2, p3 : Antione equation parameters to give vapor pressure in mmHg
        K : vapor permeability (m^2)
        h_w : outer wall newton cooling parameter (W/m^2/K)
        h_f : inner wall newton cooling parameter (W/m^2/K)
        X_cr1 : first critical solvent concentration 
        X_cr2 : second critical solvent concentration 
        z_c : cake height (m)
        D : vessel diameter (m)
        P_v : vessel pressure (Pa)
        Tf_i : initial temperature heat exchange fluid (K) # TODO: What is this?
        grad_Tf : slope change heat exchanger fluid (K/s) # TODO: What is this?
        Tf_s : set point of the heat exchanger fluid (K) # TODO: What is this?
        To : initial temperature (CK 
        Xo : initial mass fraction
        dt : integration time increment (s)
    """
    """ set initial conditions """
    t = 0                               # setting initial time (s)
    m = 0                               # evaporation rate (kg/s)
    z = 0                               # front position (m)
    X = Y/(1-Y)*r_l/r_s                 # mass fraction (based on porosity)
    T = To                              # setting starting temperature (K)
    results = [{'t':t, 'm':m, 'z':z, 'X':X, 'T':T}]

    """ calculate additional parameters """
    A = np.pi * D**2 / 4.               # cross section area of vessel (m^2)
    m_s = (1-Y)*A*z_c*r_s               # total mass of solid (kg)
    r = (P_v * mw ) / (8.314 * To)      # vaport density from ideal gas law (kg/m^3)

    """ begin constant rate phase """
    h = h_f
    while X > X_cr1:
        T_p = results[-1]['T']
        T_f = jacket_temperature(t, Tf_i, grad_Tf, Tf_s)
        [T,fopt,iter,funcalls,warnflag] = fmin(local_temp_objective_function, T, p,
                                               xtol=0.0001, ftol=0.0001, maxiter=10000, maxfun=None, 
                                               full_output=1, disp=0, retall=0)
        if warnflag != 0:
            pdb.set_trace()
        T = T[0]
        P_sat = self.vapor_pressure(T,P_parameters)
        j = self.mass_flux(P_sat,z,z_c,P_v,K,v,r)
        q = self.heat_flux(T,z,h,T_f,l_d)
        m = -1 * j * A
        X = (m/m_s)*dt + X
        self.simulation_results.append([t,m,z,X,T])
        t = t + dt
    print 'begin transition phase'
    while X > X_cr2:
        h = self.Newton_cooling_parameter(X,h_w,h_f,X_cr1,X_cr2)
        T_p = self.simulation_results[-1][4]
        T_f = self.calculate_jacket_temperature(Tf_i,grad_Tf,Tf_s,t)
        p = (P_parameters,z,P_v,K,v,r,h,T_f,H_vap,z_c,l_d,cp_s,r_s,cp_l,r_l,X,T_p,dt)
        #data = self.regression_diagnostics_inspection(numpy.linspace(T-20,T+20,100),*p)
        [T,fopt,iter,funcalls,warnflag] = fmin(self.local_temp_objective_function,T,p,
                                               xtol=0.0001,ftol=0.0001,maxiter=10000,maxfun=None,
                                               full_output=1,disp=0,retall=0)
        if warnflag != 0:
            pdb.set_trace()
        T = T[0]
        P_sat = self.vapor_pressure(T,P_parameters)
        j = self.mass_flux(P_sat,z,z_c,P_v,K,v,r)
        q = self.heat_flux(T,z,h,T_f,l_d)
        m = -1 * j * A
        X = (m/m_s)*dt + X
        self.simulation_results.append([t,m,z,X,T])
        t = t + dt
    print 'begin falling rate phase'
    while X > 0:
        T_p = self.simulation_results[-1][4]
        T_f = self.calculate_jacket_temperature(Tf_i,grad_Tf,Tf_s,t)
        p = (P_parameters,z,P_v,K,v,r,h,T_f,H_vap,z_c,l_d,cp_s,r_s,cp_l,r_l,X,T_p,dt)
        #data = self.regression_diagnostics_inspection(numpy.linspace(T-20,T+20,100),*p)
        [T,fopt,iter,funcalls,warnflag] = fmin(self.local_temp_objective_function,T,p,
                                               xtol=0.0001,ftol=0.0001,maxiter=10000,maxfun=None,
                                               full_output=1,disp=0,retall=0)
        if warnflag != 0:
            pdb.set_trace()
        T = T[0]
        P_sat = self.vapor_pressure(T,P_parameters)
        j = self.mass_flux(P_sat,z,z_c,P_v,K,v,r)
        q = self.heat_flux(T,z,h,T_f,l_d)
        m = -1 * j * A
        X = (m/m_s)*dt + X
        z =  j /(X_cr2*(1-Y)*r_s) * dt + z
        X = X_cr2 * (z_c - z) / z_c
        self.simulation_results.append([t,m,z,X,T])
        t = t + dt
    print 'evaporation completed'
    print 'beginning transient heating'
    self.transient_heat_conduction()
    print 'simulation complete'
    return numpy.array(self.simulation_results).transpose()

def transient_heat_conduction(self):
    """
    Calculate the transient heat up of the dry layer
    """
    l_d = self.l_d 
    cp_s = self.cp_s 
    r_s = self.r_s 
    z_c = self.z_c
    Tf_s = self.Tf_s - 273
    dt = self.dt
    dz = z_c/100.
    r = l_d*dt/(r_s*cp_s*dz)
    T = numpy.zeros((2,100))
    T[0] = numpy.linspace(Tf_s,self.simulation_results[-1][4],100)
    T[:,0] = Tf_s
    previous_results = self.simulation_results[-1][1:4]
    tm = self.simulation_results[-1][0]
    for t in range(1,10000):
        for z in range(1,99):
            T[1,z] = r*T[1-1,z-1] + (1-2*r)*T[1-1,z] + r*T[1-1,z+1]
        T[1,99] = 2*r*T[1-1,z-1] + (1-2*r)*T[1-1,z]
        T[0] = T[1]
        tm = tm + dt/60.
        if numpy.mod(t,100)==0:
            self.simulation_results.append([tm]+previous_results+[T[1,99]])

def make_plots(self,data=None):
    if data == None:
        data = numpy.array(self.simulation_results)
    fig = pyplot.figure()
    ax = fig.add_subplot(111)
    ax.plot(data[:,0],numpy.abs(data[:,3]),color='blue',linestyle='-')
    ax.set_xlabel(r'time (s)')
    ax.set_ylabel('mass fraction (kg solvent/kg dry solid)')
    ax2 = ax.twinx()
    ax2.plot(data[:,0],data[:,4],color='green',linestyle='-.')
    ax2.set_ylabel(r'drying front temperature (C)')
    fig = pyplot.figure()
    ax = fig.add_subplot(111)
    ax.plot(data[:,3],numpy.abs(data[:,1]),color='blue',linestyle='-')
    ax.set_xlabel('mass fraction (kg solvent/kg dry solid)')
    ax.set_ylabel(r'drying rate (kg solvent/s')
    pyplot.draw()

if __name__=='__main__':
    pass

